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We carried out a multiwavelength cross-correlation analysis of a sample of 16 blazars detected by 
Fermi/LAT. The purpose is to investigate if there exists correlations between the distinct bands we 
analyze in this work. We searched for cross-correlated delays between emissions in optical, near-infrared 
and 7-ray bands for 16 blazars detected by Fermi-LAT, using three methods previously discussed in 
the literature: the interpolated cross-correlation function, the discrete cross-correlation function and 
the Z-transformed discrete cross-correlation function. Our results confirm the expectation that in all 
our sample the four NIR bands vary simultaneously. For three objects of our sample (3C 273, Mrk 501, 
and PMN J0808-0751), no correlation is found between any of the bands available for this study. For 
the remaining thirteen, a correlation was found between the V band and the NIR bands, indicating 
that in most of them the V band and the NIR bands vary simultaneously. For 4 objects (3C 454.3, 
PKS 0235+164, PKS 1510-089, and PKS 2155-304) a delayed correlation between the 7-ray emission 
and the NIR+V bands was found. 
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1. INTRODUCTION 

Blazars are a subclass of Active Galac- 
tic Nuclei (AGN) characterized by a promi- 
nent jet whose angle relative to the line of 
sight is very sma l l, wi thin a few degrees 



Blandford fc Rees 1978 



In 



many cases, 



these jets can be very luminous at wavelengths 
from radio to 7-rays. Due to the proximity 
of the jet axis to the line of sight, the emis- 
sion from the jet is relativistically beamed, 
and hence amplified by an order of magnitude 
or more for some blazars. The observed time 
scales, in some cases, are also shorter than 
those in the rest frame of the jet. Most of 
the observed emission from radio to optical 
(sometimes UV) is due to synch rotron radia- 
tion from the jet |Marscherl 1998 |. X-rays and 
7-rays may be produced via Inverse Comp- 
ton scattering by the same energetic electrons 
radiating synchrotron em ission (the so-called 



lepto nic models, e.g. iBottcher fc Dermer 
20071 ) ; or may be due to synchrotron radiation 
by protons co-accelerated with the electrons 
in the jet, interactions of these highly rela- 
tivistic protons with external radiation fields, 
or proton-in duced particle casc ades (hadronic 



models, e.g. iMiicke et al.l 120031 1 . 



Time variability at multiple wavelengths is 
a defining characteristic of these objects, and 
has been used to probe the location and phys- 



ical processes related t o the emission at ver y 
fine resolutions (e.g. IMarscher et al.l 120081 ). 
Most of the observed blazars have the peak 
of their emission in the 7-ray part of the 
spectrum. However, until recently, long-term 
and well-sampled 7-ray light curves of blazars 
were not available. This has changed with 
the near-continuous monitoring activity of the 
Large Area Telescope (LAT) instrument on 
board the Fermi Gamma-Ray Space Tele- 
scope, launched in 2008, which provides the 
opportunity to study the variable SEDs of a 
large sample of blazars with truly simultane- 
ous multi-frequency data. 

In this contribution we present the results 
of a multiwavelength variability study on a 
sample of blazars observed by Fermi/LAT. 



2. THE SAMPLE AND OBSERVATIONAL 
DATA 

A sample of 16 blazars detected by 
Fermi/LAT was analyzed: 9 Optically Violent 
Variable quasars (OVV) and 7 BL Lacertae 
type objects (BLL) (see Table HI). 

We have light curves in Near InfraRed 
(NIR) J, H, K and Ks bands; and optical 
(V band) monitoring of Fermi/LAT Blazars. 
The time period covered by the observations 
is from the late 2007 up to mid-2012. The 
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Table I The Sample 



Object 


NED Classification Redshift 


V 


3C 66A 


BLL 


0.444 


14.4 


BL Lac 


BLL 


0.069 


14.7 


Mrk 421 


BLL 


0.030 


13.0 


Mrk 501 


BLL 


0.034 


13.9 


PKS 0716+714 


BLL 


0.300 


13.6 


PKS 2155-304 


BLL 


0.116 


13.5 


W Comae 


BLL 


0.102 


15.2 


3C 273 


OVV 


0.158 


12.8 


3C 279 


OVV 


0.536 


16.5 


3C 345 


OVV 


0.593 


17.2 


3C 454.3 


OVV 


0.859 


15.6 


PKS 0235+164 


OVV 


0.940 


18.0 


PKS 1510-089 


OVV 


0.360 


16.5 


PKS 1633+382 


OVV 


1.814 


17.4 


PMN J0808-0751 


OVV 


1.837 


R~16.7 


QSO B0133+476 


OVV 


0.859 


17.0 



NIR data are from the Guillermo Haro As- 
trophysical Observatory (OAGH from Span- 
ish) using the Cananea Near-Infrared Camera 
(CANICA) and f rom from the Yale Fermi/ 
SMARTS project |Bonning et al.ll2012j |: while 
the optical data (V Band) are from four dif- 
ferent sources: most of the data is from the 
Ground-based Observational Support of the 
Fermi Gamma-Ray Space Telescope a t the 
University of Arizona (Smith et al.ll2009 | using 
the SPOL CCD Imaging/Spectropolarimeter 
and from the Yale Fermi/SMARTS project. 
Extra data for 3C 454.3 was taken from 
Raiteri et all |201l| . the WEBT Project; and 
V band data for QSO B0133 +476 wa s taken 
from the MIS AO Project [Yoshidal 12012 



specifically photometry taken by Miguel Ro- 
driguez Marco with the SRO50 AAVSONet 
Robotic Telescope. 

3. CROSS CORRELATION ANALYSIS 

A Cross-Correlation analysis of the light 
curves was carried out with the pur- 
pose of quantify lags between different 
kinds of emission, using three different 
methodologies of the cross correlation func- 
tion: The Interpolated Cross Correlation 



Function (ICCF, iGaskell fc Sparkd Il986h . 
The Discrete Cross Correlation Function 



(DCCF, lEdelson fc Krolikl Il988h and the Z- 
Transformed Discrete Cross Correlation Func- 



tion fZDCF. lAlexanderl 119971 ). 



3.1. The Interpolated Cross Correlation 
Function 



The ICCF |Gaskell fc Sparkd Il986| ] consists 
in interpolate the actual light curves, so to get 
the same number of points and equal sampling 
in both light curves. Once we have vectors of 
identical size, covering a time interval T for 
both light curves, we must define a vector of 
lags which can run from — T to +T (in units 
of the spacing of the interpolated points); this 
vectors of lags represent movements that will 
have the time axis of one of the light curves. 
When we move the light curves in time, we 
compute the correlation coefficient only be- 
tween the segments of the light curves that 
intersect; this is done to avoid problems for 
the non-stationarity of the light curves (see 
White fe Peterson! 1 1994T ). Note that if we dis- 
place the time axis of one of the light curves by 
a lag +T or — T, we will be left with just one 
point in each light curve to compute the cor- 
relation coefficient, which is not statistically 
significant; therefore we make the vector of 
lags to run from — 0.8T to +0.8T so we are 
left with at least 20% of the light curve to 
compute a correlation coefficient. 

Once we have our light curves x and y, 
and a vector of lags L, we can calculate the 
Cross-Correlation Func tion P T JL) by using 



the Eq. [T] as defined by iFullerl 1976I ] : 



N-\L\-1 

( X k+\L\ - x)(yk - y) 



k = 



Pxy(L) = ■ 



N-i "| nv-i 

Y ( x k ~ x ) 2 ^ k ~ 

fe=0 J Lk=0 

N-L—l 

Y (a;* -x)(y k+L -y) 
fc=0 



For L < 



For L > 



N-l 
k=0 



N-l 
fc=0 



(1) 



eConf C121028 



4 th Fermi Symposium : Monterey, CA : 28 Oct- 2 Nov 2012 



3 



10 

"« 8 

" 6 

a 

o 

I 4 
o 

x 2 


16 



1 



12 
g 
4 
80 

60 

' 40 

20 



1 I r " 

Fermi - y 0.2 - 300 GeV 



n I 
3C 454.3 




Steward + WEBT + SMARTS - V 



•i tit ! 



OAGH + SMARTS - J 



. i 



4500 5000 5500 

JD - 2450000 



6000 



S 6 



5 2 


5 

4 

>. 3 
1 

2 
1 


16 

12 



TTn 1 i 1 r 

. Fermi - y 0.2 - 300 GeV 3C 279 



». 



Steward + SMARTS - V 



» . 



OAGH + SMARTS - J 



7 • 



4500 5000 5500 6000 

JD - 2450000 



Figure 1: Light Curves from 3C 454.3 and 3C 279 in 7-rays, V band and J band. Light brown squares in the 3C 454.3 
7-rays light curve are upper limits. We also have emission line and spect ral continuum light curves from the Steward 
Observatory, which are being complemented with data from OAGH (see iPatino- Alvarez et al.ll2013T l. 



where x is the mean of the points in the x 
curve that are involved in the calculation of 
the correlation coefficient, and similarly for y. 
The means are recalculated for each lag. 

3.2. The Discrete Cross Correlation 
Function 

Unlike the ICCF, the method proposed by 
Edelson fc Kroliki 1988| does not assume any 
understanding of the real light curve behavior; 
it uses only real data points which are sepa- 
rated by the time r. The correlation function 
itself is binned on time intervals St, so that 
the value of the DCF at r is the average over 
the interval r — St/2 to r + St/2. 

If you have two vectors a and b, then you 
can form pairs (a*, fy), each one of these is 
associated with the pairwise lag At^ = tj —U. 
Now, for each of these pairs we can compute 
the unbinned discrete correlation functions: 



UDCR, 



(cii - a) (bj - b) 



(2) 



where a and a a are the means and standard 
deviations of the points in each bin, and in 
similar manner for series b. 



Originally lEdelson fc Kroliki [1988| used the 
entire series to calculate the means and the 
standard deviations, which is only applicable 
when we have stationary light curves, how- 
ever, the AGN light curves are not station- 
ary (i.e. the mean and standard deviation 
change with time 



199 



e.g. 



White fc Peterson 



Also, they suggested that in order to 
preserve the proper normalization it is nec- 
essary to replace the a a &b in the Eq. [2] with 



Keg - e 2 a )(a 2 - el)} 1 ' 2 . IWhite fc Peterson! 
|1994| . however, argue that doing this greatly 
complicates any direct comparison with the 
interpolation results, which are unweighted, 
and therefore this makes an interpretation of 
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the DCF amplitude less straightforward. 

Next, if we have M pairs for which r — 
St/2 < Aty < t + St/2, then the discrete 
correlation function of r is: 



DCF(r) 



1 

M 



•rSt/2 



UDCF^Atij) (3) 

T-St/2 



[Note that the DCF(r) is not defined for a 
bin with no points] 

Unlike the interpolation method, for which 
the errors are difficult to define, we can de- 
fine almost directly a standard error for the 
DCF. If each of the individual UDCF^ within 
a single bin were totally uncorrelated, then the 
standard error in the determination of their 
mean would be: 



cdcf(t ) = 



•JM{M - 



Y^UDCFij ~DCF(t) 



1/2 



(4) 



3.3. The Z-Transformed Discrete Cross 
Correlation Function 

The Z-Transformed Discrete Correlation 



Function by lAlexanderl [1997j is an alterna- 
tive method for estimating the CCF of sparse, 
unevenly sampled light curves. The ZDCF 
corrects several biases of the discrete corre- 



lation function method of lEdelson fc Krolik 
19881 ] by using equal popul ation binning and 



(e.g. 



Kendall fe Stuart 



Fisher's z- transform 
1969 . 19731 and references therein). 

If we have n pairs in a given time-lag bin, 
the CCF(t) is estimated by the correlation 
coefficient: 



Y, n M-a){h-b)/{n-l) 



SaSb 



(5) 



where a, b are the bin averages, and s a , Sb 
are the standard deviations of the points in a 
bin. 



The sampling distribution of r is highly 
skewed and far from normal, therefore esti- 
mating its sampling error by the simple vari- 
ance s r can be very inaccurate. 

If a and b are drawn from bivariate nor- 
mal distributions, it is possible to transform 
r into an approximately normally distributed 
random variable, Fisher's z. Defining: 



z — — m I I , p = tann z , (, = — m 



2 1-r 



2 I 1 -p 



(6) 



This yields to the mean of z being approxi- 
mately equal to: 



2(n - 1) 



5 + p 2 ll + 2p 2 +3p 4 
4(n-l) 8(n-l) 2 



(7) 



and the variance of z is approximately equal 



to: 



n - 1 



4 - p 2 22 - 6p 2 - 3p 4 



2(n - 1) 



6(n- l) 2 



(8) 



Transforming to r again, the interval corre- 
sponding to the normal ±lcr error interval can 
be determined by: 



Sta 



tanh (z(r) ± s z (r)) — p| (9) 



The binning met hod for the ZDCF is dif - 



1988 



ferent from that of lEdelson fc Krolik 
in which the binning of the pairs is done in 
a fixed time interval, i. e. every bin has the 
same length in time, and the resulting DCF 
of that bin is the average of the UDCF of all 
pairs who fall into the bin; while, the ZDCF 
binning is for a fixed population, i. e. ev- 
ery bin has the same number of pairs, at least 
n min = 11 pairs, which is the minimum num- 
ber for a meanin gful statistical interpretation 



Alexander 1997 



Also, in each bin, the interdependent pairs 
are discarded. This means that in a bin we 
can not have two pairs which use the same 
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5 



di or bj element. Therefore, light curves with 
less than 12 points cannot be analyzed by this 
method. 

We made a few modifications to the original 
method, by changing how the accommodation 
of pairs is done, we do not start from the me- 
dian. Since many of the pairs found in the 
median of the A^- distribution use some of 
the same points of the light curves, a lot of 
pairs are discarded and not used in the actual 
analysis; therefore we start the accommoda- 
tion of the pairs from the extreme negative 
part of the CCF (i.e. the pair with the small- 
est Aty). Also, we require that all bins have 
the same number of pairs, unlike the original 
method where it was possible to obtain bins 
with different number of pairs; which may lead 
to different significance in each bin. 

Examples of the CCF curves obtained with 
the three different statistical methods are 
shown in Fig. [2] 

4. RESULTS 

For the objects in our sample we have light 
curves in 7-rays, V band, J, H, K and Ks 
bands: 7-rays and NIR light curves for all of 
them, one object doesn't have a V-band light 
curve. The lags obtained between the differ- 
ent bands in all objects can be found in Ta- 
ble Ull The main conclusions about the sample 
can be summarized as follows: 

1. We found that in all objects the four NIR 
bands (J, H, K and Ks) vary simultane- 
ously. 

2. No significant difference was found in the 
results obtained for the OVV, and the 
results obtained for the BLL. 

In the specific case of the source 3C 279, 
the preliminary analysis of different parts of 
the available light curves, showed two differ- 
ent significant peaks in the cross-correlation 
functions obtained with the three methods; 
this could mean that we have two different 
delays between the 7-rays emission and the 
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Figure 2: Examples of the CCF curves obtained with the 
three different methods. 

optical/NIR emissions: One near zero days, 
and the other in ~ 250 days. 

Our results that in most sources the NIR 
and optical emissions are well correlated is in 
agreement with the results obtained by other 
authors. Our results also confirm the corre- 
lations found for PKS 1510-089 and 3C 454.3 
between 7-rays and the optical/NIR emissions 
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e.g. iBonning et al.l 120121 . iMarscher et al. 



20101 ) . 

The multiwavelength analysis will be com- 
plemented with millimeter and X-rays data. 
More details about the analysis will be in a 
forthcoming paper. 

Table II General results of the Cross-Correlation analysis. 



Object 


7-rays vs. V 


7-rays vs. NIR V vs. NIR 


3C 66A 


N/C 


N/C 


11 7+ 7 - 1 

' -7.0 


BL Lac 


N/C 


N/C 


-79 ± 13 


Mrk 421 


N/C 


N/C 


3 4+ 69 


Mrk 501 


N/C 


N/C 


N/C 


PKS 0716+714 


N/C 


N/C 


-1 ± 13* 


PKS 2155-304 


3.5 + 2.3 


3.7 + 2.4 


-o.7±l:f 


W Comae 


N/C 


N/C 


r\ -i +8.3 

^■ L -8.2 


3C 273 


N/C 


N/C 


N/C 


3C 279 


** 


** 


0+ 31 


3C 345 


N/C 


N/C 


J-' -12 


3C 454.3 


1.1 ± 1.8 


-0.4+1.8 


0.0 ± 1.2 


PKS 0235+164 


221J 1 


18+14 


1 ± 12 


PKS 1510-089 


7.3t«* 


22 ± 10 


9+11 


PKS 1633+382 


N/C 


N/C 


O -i +8.5 
z - i -4.5 


PMN J0808-0751 


N/V 


N/C 


N/V 


QSO B0133+476 


N/C 


N/C 


-4 5+ 9 5 

9.1 



These results are uncertain, due to problems during 
analysis. 

Probable two delays, deeper study needed. 
N/C: No correlation found. N/V: No V band. 



Acknowledgments 

This work was supported by CONACyT re- 
search grant 151494 (Mexico). V.P.-A. ac- 
knowledges support from the CONACyT pro- 
gram for PhD studies. 



References 

Alexander, T. 1997, Astronomical Time Se- 
ries, 218, 163 

Blandford, R. D., & Rees, M. J. 1978, in 
Pittsburgh Conference on BL Lac Objects 



Proceedings. (A79-30026 11-90), ed. A. M. 
Wolfe (Pittsburgh, Pa., University of Pitts- 
burgh), 328-341 

Bonning, E. W., Urry, C. M., Bailyn, C, et al. 
2012, ApJ, 756, 13 

Bottcher, M., & Dermer, C. D. 2007, Ap&SS, 
307, 233 

Edelson, R. A., & Krolik, J. H. 1988, ApJ, 
333, 646 

Fuller, W. A. 1976, Introduction to Statistical 
Time Series (John Wiley & Sons Inc.) 

Gaskell, C. M., & Sparke, L. S. 1986, ApJ, 
305, 175 

Kendall, M. C, & Stuart, A. 1969, The Ad- 
vanced Theory of Statistics, 3rd edn., Vol. I 
(Griffin, London) 

— . 1973, The Advanced Theory of Statistics, 
3rd edn., Vol. II (Griffin, London) 

Marscher, A. P. 1998, in in ASP Conf. Ser., 
Vol. 144, Radio Emission from Galactic and 
Extragalactic Compact Sources, ed. J. A. 
Zensus, J. M. Wrobel, & G. B. Taylor (San 
Francisco, CA: ASP), 25 

Marscher, A. P., Jorstad, S. G., D'Arcangelo, 
F. D., et al. 2008, Nature, 452, 966 

Marscher, A. P., Jorstad, S. G., Larionov, 
V. M., et al. 2010, ApJ, 710, L126 

Miicke, A., Protheroe, R. J., Engel, R., 
Rachen, J. P., & Stanev, T. 2003, APh, 18, 
593 

Patino-Alvarez, V., Chavushyan, V., Leon- 
Tavares, J., et al. 2013, in 4th. Fermi 
Symposium, eConf Proceedings C121028. 
(arXiv: 1303. 1893) 

Raiteri, C. M., Villata, M., Aller, M. F., et al. 
2011, A&A, 534 

Smith, P., Montiel, E., Rightley, S., et al. 
2009, in 2nd. Fermi Symposium, eConf Pro- 
ceedings C091122. (arXiv:0912.3621) 

White, R. J., & Peterson, B. M. 1994, Publ. 
Astron. Soc. Pac, 106, 879 

Yoshida, S. 2012, the MISAO project (online 
data) 



eConf C121028 



